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Abstract 

We present Tethered Monte Carlo, a simple, general purpose method of computing the efTec- 
tive potential of the order parameter (Helmholtz free energy). This formalism is based on a 
new statistical ensemble, closely related to the micromagnetic one, but with an extended con- 
figuration space (through Creutz-like demons). Canonical averages for arbitrary values of the 
external magnetic field are computed without additional simulations. The method is put to work 
in the two dimensional Ising model, where the existence of exact results enables us to perform 
high precision checks. A rather peculiar feature of our implementation, which employs a local 
Metropolis algorithm, is the total absence, within errors, of critical slowing down for magnetic 
observables. Indeed, high accuracy results are presented for lattices as large as L = 1024. 
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1. Introduction 

There exists a very profound relation between quantum field theory and statistical 
mechanics, through the theory of critical phenomena [1-3]. In the examination of these 
phenomena Monte Carlo simulations [4] are an indispensable tool, because of their ex- 
ceedingly broad range of applicability. Indeed, Monte Carlo simulations not only succeed 
where an analytical treatment would be impossible or impractical, but they are many 
times the only numerical method capable of tackling the problem at hand. 
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There are some difficulties, though, among which we can mention critical slowing 
down [2,5]. For traditional Monte Carlo formalisms the correlation times (roughly, the 
number of intermediate steps so that two measurements can be considered independent) 
grow as at the critical point (with 2; « 2). A great step towards the solution of this 
problem was taken in the late 1980s, with the development of the first cluster methods [6], 
capable of achieving z < 1. This prompted a large amount of work on more sophisticated 
cluster algorithms, which continues today [7]. 

Unfortunately, cluster methods are highly specific and wc do not know how to efficiently 
implement them for physical problems as important as lattice gauge theories [8], with 
or without dynamic fermions; structural glasses [9]; spin glasses [10]; protein folding [11] 
and a long etcetera. Even in their favourite playground, ferromagnetic systems, chistcr 
methods lose most of their power in the presence of a magnetic field. The simplest example 
is the D — 2 Ising model, whose exact solution without a magnetic field is known since 
1944 [12], but whose behaviour with a magnetic perturbation is still an active research 
topic [13,14]. It is interesting to notice that the current numerical methods of choice rely 
on transfer matrix techniques [15] and not on Monte Carlo simulations. 

Here we present the Tethered Monte Carlo (TMC) method, a completely unspecific 
strategy, appliable to many problems with or without an external field. The main goal 
of this approach is constructing the effective potential of the order parameter (perhaps 
more commonly named Helmholtz free energy in a statistical mechanics context). 

In order to define this Monte Carlo method, we shall introduce a new statistical ensem- 
ble, where the magnetic field and the order parameter exchange their roles with respect 
to the canonical ensemble. ^ This new framework has some interesting features of its 
own. For example, when working with a ferromagnetic system it provides a very clean 
definition of the broken symmetry phase for a finite lattice. 

In the ferromagnetic setting the tethered ensemble is related to the micromagnetical 
one, where both the temperature /3 and the magnetisation density m are kept fixed. The 
difference is that here we couple m to a Gaussian 'magnetostat'. This yields a new variable 
m which, unlike m, is continuous even for finite lattices. The magnetic field is obtained 
from a fluctuation-dissipation formalism. We can then work at constant m, where the real 
magnetisation is almost fixed, but has some leeway (hence the name tethered). We then 
combine the mean values at constant m of the magnetic field to construct the effective 
potential i7jv(m), whose exponential gives the canonical probability density function 
(pdf) of m. Using this pdf and the mean values of the different physical observables as a 
function of m we can recover the canonical results. 

The tethered ensemble is not only a fancy theoretical construct, but also the basis for 
a practical simiilation method. It is straightforward to implement it with, for example, 
the well known local Metropolis [16] update algorithm. We can thus run simulations 
at constant m and then combine them to obtain the canonical averages with very high 
precision. It would be also possible to search for a cluster algorithm compatible with the 
tethered formalism, but we shall not investigate this issue here. 

One surprising feature of this Metropolis implementation is that we can find no traces 
of critical slowing down, within our errors, for all functions of m. Other quantities, such 
as the energy, exhibit the 2 « 2 behaviour typical of local algorithms. Another interesting 



^ In the canonical ensemble the magnetic field is an external parameter and the order parameter an 
observable, while their tethered equivalents have the reversed roles. 
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point about this method is that, for a given temperature, f2N{'ni) has all the information 
about the system so we can, for example, extract the canonical results at any value of 
the applied magnetic field h without recomputing the tethered mean values (i.e., without 
any new simulations, see Section 7 below). 

The TMC method is an extension of the strategy introduced in [17]: there the config- 
uration space was extended to work in a microcanonical ensemble, with entropy as the 
main physical variable. The motivation in [17] was handling first order transitions with- 
out the need for tunnelling between two coexisting phases. The microcanonical method 
was first applied in [17] to a pure system and was later, in [18], instrumental in simulat- 
ing a phase transition which remained first order in the presence of quenched disorder. 
Both the TMC and microcanonical methods have deep links to Crcutz's microcanonical 
demon [19]. The main differences are: (i) we have as many demons as spins, (ii) our 
demons are continuous variables and (iii) we explicitly integrate out the demons, finding 
a tractable effective Hamiltonian. 

The effective potential can also be computed using multicanonical methods [20] , some- 
times named multimagnctical [21], or with the Wang-Landau algorithm [22]. A major 
difference is that TMC does not require a random walk in magnetisation space (as with 
multimagnctical methods) or in the energy-magnetisation plane (as in Wang-Landau). 
Indeed, we have worked with as many as 10® spins, while 10"^ spins is a typical limit for 
Wang-Landau computations [23]. On the other hand, standard micromagnetical meth- 
ods [24] do not render the effective potential. 

In this paper we give a detailed exposition of the TMC method and demonstrate it 
in the standard benchmark of the two dimensional Ising model. Our motivation for this 
choice is twofold. On the one hand, since Onsager's solution [12] many other exact results 
have been obtained (for a review see [25] and references therein), which will help us check 
that our answers are correct. This model is also the ideal scenario for cluster methods, 
so even for those observables whose exact value is unknown we can check our results to a 
high degree of accuracy against a cluster simulation. On the other hand, the Ising model 
is sufficiently well known and simple that we may concentrate on examining the details of 
the method with a minimum of nonessential discussion. We purpose to show that TMC 
is capable of rendering very precise results and try to convince the reader that it will still 
be efficient for harder problems. 

The rest of the paper is organised as follows: 

• In Section 2 we describe our statistical ensemble and its relationship with the standard 
canonical ensemble and properly define the effective potential. In Section 3 we explain 
in detail how to set up a simulation using the TMC method. We then examine our 
own simulations for the Ising model. In Section 4 we show that our algorithm presents 
no measurable traces of critical slowing down for the magnetic field. 

• We have carried out simulations at the critical point and in both the ferromagnetic 
and paramagnetic phases, with and without a magnetic field. Our results at the critical 
temperature and zero magnetic field are presented in 5 and compared both with the 
exact results at finite L given by [26] and with high precision Swendsen-Wang simula- 
tions (Section 5.2). In Section 5.3 we shall compute very accurately the magnetisation 
critical exponent in a fairly unusual way, made easy by TMC. In Section 6 we check 
the performance of the method in the scaling paramagnetic region. Our chosen test 
has been the computation of the first renormalised coupling constants. 

• The magnetic field h is introduced in Section 7, where we obtain several observables 
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as a function of h and check our results with FSS arguments and by recomputing the 
nonlinear susceptibilities with a finite differences formula. Section 8 demonstrates the 
enhanced effectivity of the method in the ferromagnetic regime, where we work with 
and without an external field. 
• Finally, in Section 9 we present our conclusions and outlook. We give some technical 
details of our numerical implementation in an Appendix. 



2. The Statistical Ensemble 



2.1. The model and the canonical observables 



We shall work with the D = 2 Ising model, defined by the following partition function, 
Z= 5^ e''^<-«> (T^=±l, (1) 

where h is the applied magnetic field, x = {xi,X2) and {x,y) stands for lattice nearest 
neighbours. We shall always consider square lattices with N = spins and periodic 
boundary conditions. The infinite volume model has a second order phase transition at 
a critical (inverse) temperature (3c given by 

/J^ = i^iii^l^ = 0.440 686 793 509 771... (2) 
The simplest observables are the energy and magnetisation, ^ 

U = Nu = -Y, a^ay, (3) 
{x,y) 

M = Nm = J2<^x. (4) 

X 

In the canonical ensemble we are concerned with thermal averages, which we shall denote 
by ( • 

U0 = NU0 = {U)0, (5) 
= Nmp = {M)0. (6) 

The specific heat and magnetic susceptibility are 

C = N[{u^),-{u)ll (7) 
X2 = N[{m^}^-{m)l]. (8) 

The latter quantity can be seen as the zero momentum component of the two point 
correlation function: 

G2{k) = j^Y.(^xao)0 e"=-. (9) 



^ We shall use sans-serif italics for random variables (i.e., functions of the spins) and serif italics for 
real numbers (e.g. expectation values or arguments of the probability density functions). This will help 
emphasise which quantities are kept fixed and which can change. We shall also use lowercase letters for 
intensive quantities and uppercase letters for extensive quantities. 
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If we consider the asymptotic behaviour of the propagator in position space we arrive at 
the concept of correlation length, 



Sexp 



lim 

x\ — *cx 



(10) 



This quantity is not easily measurable in our finite lattices, so we would like to have 
a better statistically behaved definition that could also be interpreted as a correlation 
length. In momentum space and in the limit ^exp|fc| 0, the propagator is well described 
by the free field form [1,2] 

{A is a constant). If we combine this formula at zero momentum and at the smallest 
nonzero momentum fei we obtain 

1/2 

(12) 



6 = 



1 



2sin(7r/L) 



G2(0) 



G2(fel) 



- 1 



with G2{ki) averaged over ki = (27r/L,0) , {0,2tt/L). This definition [27] has proven to 
be extremely useful in Finite Size Scaling studies [28,29]. 

Equation (12) does not work in the broken symmetry phase or with an applied magnetic 
field, because then G2(fe) becomes singular at A; = 0. Wc can still use (11) for fe 7^ and 
consider a second definition of the correlation length, now combining the two smallest 
nonzero momenta, 



1 



G2(fcl)-G2(fe2) 



-,1/2 



2sin(Vi) [2G2(fc2)-G2(fci)J ' ^^^^ 

with G2(fc2) averaged over k2 = {2^ / L , ±2n / L) and G2{k\) as in definition (12). The 
two definitions, ^1 and ^2, coincide for /3 < /3c, but only in thermodynamic limit. 
Other observables are the Binder ratio 



B = 



'/3 



(M 



2n2 

//3 



and the 2n point correlation functions at zero momentum, 

1 52"log(Z) 



X2n 



N {dh) 



2n 



(14) 



(15) 



Notice that these are just the cumulants of the magnetisation. In the low temperature 
phase we should also consider odd derivatives. 



2.2. The tethered ensemble 

Let us consider the Ising model without an external field {h = 0). We can define a 
probability density function (pdf) for M, which will be a sum of A'^ + 1 Dirac deltas, 

Pi{m) = I E eM-PU]s(m -Y,^i/N\ , Z=Y, exp[-/3(;]. (16) 
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In the thermodynamic hmit pi{m) is non vanishing for all m in [—1, 1]. Wc would like to 
have a new quantity that would be continuous even for finite lattices. As a first step, we 
extend our configuration space with N Gaussian demons: ^ 



n drn E '^^P -f^U - E ^'/2 , R = Nr = Y^ vV^- (17) 

i=l {a^} L ^ J ^ 

P2{r) = ^ E<^^pf-/5f^-E^H^("-E^'/(2^))- (18) 

Notice that in the canonical formalism the demons are a thermal bath decoupled from 
the spins because the spin contribution in (18) cancels out. By virtue of the Central Limit 
Theorem, the pdf for r approaches a Gaussian of mean 1/2 and variance {2N)~^ in the 
large N limit. 

Now we introduce M = Nrh = M + R. The new variable m is continuous and its pdf 
is simply the convolution of those of m and r (as these are statistically independent): 



p{rh) = j dm j dr pi{m)p2{r)S{m — m — r). (19) 

Notice that M > M and that p(m) is essentially a smoothed version of pi{m — 1/2). 
Finally, we introduce what will be our basic physical quantity, the effective potential 
QN{m,l3), 

p{m) = exp[A^i7Ar(m,/3)]. (20) 

The effective potential has all the information about the system at inverse temperature P, 
including what would happen if it were immersed in an external magnetic field. 

Our next step is constructing the statistical ensemble and its relationship with the 
canonical one. In order to do this we represent f2]s[{m,P) as a functional integral and 
integrate out the demons, 

N 



^Na.iA,p) = i f[d,7. E e-''''-^.'''/'<5(m - m-E»?'/(2iV) 
= I n E (m-m-Y: vVim) 



The condition m > m is explicitly enforced by the Heaviside step function 6. Notice 
as well that wc have constructed the effective potential from the starting point of a 
canonical ensemble. It would be elementary to retrace our steps for a microcanonical fi^ 
(we would just have to change the exponential of the energy in the previous equations 
to the appropriate microcanonical weight, see [17]). 



^ This is not the only possible choice. In fact, we have also experimented with Poissonian demons, 
better suited to a possible future implementation of this method in dedicated computers with specific 
hardware [30]. The results (both in simulation time and accuracy of the final values) were virtually 
identical. 
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We want to develop a statistical ensemble based on the effective potential (i.e. to define 
average values). To do this we differentiate Qn, 

a/?^(m, /?) _ (-1 + 2lv^^) m, N- {<tJ) ^^^^ 



dm Y.{a^}^{l^,rn,N;{a^}) 



where 



u{[3, m, N- {cr^}) = Q-0U+M-Nrh^y^ _ ^)(JV-2)/2 _ (23) 

Equations (22) and (23) suggest a new ensemble where the probability of a given config- 
uration {(Jx} is proportional to Lo{P,rh, N; {dn}). Therefore, we can define the tethered 
mean value { ■ )m,i3 for an arbitrary observable by 

E{a.} 0(m; {c7^})w(/3, m, TV; {(Ta,}) 

Now we define the tethered magnetic field as 

h(rh;{a^}) = -l+ ^P~^ (25) 
M-M 

and notice from (22) and (24) that 

{hU,^ = • (26) 

The duality between the roles of h and m in the canonical and tethered ensembles is 
best illustrated from the tethered fluctuation-dissipation formula: 



d{OU,0 _/dO\ 
dm \dm/.g 



{Oh)rh.,l3 - {0)rh,(i{h)rh,(i ■ (27) 



The simulation strategy is then clear: wc compute the tethered averages of h and 
whichever observables we need for a reasonable number of values of m (we shall discuss 
what we mean by 'reasonable' in the next section). The effective potential cannot be 
measured directly, but wc do have its derivative {h)m.p- Integrating /3 and requiring 
that the probability p(m) be normalised we can obtain fiN{m,P) unambiguously. 

Once we have the effective potential, we can recover the canonical averages with the 
following formula: 

{0)p = j dm (0)a,/3 e^^^*''^) . (28) 

Thus far, we have defined the ensemble in the absence of an external magnetic field ft-, 
but we can include it very easily. Indeed, we just have to notice that an applied field 
introduces a shift in the origin of h. The computation of canonical averages with an 
external field is then straightforward, using the same tethered mean values we had for 
/i = 0: 

/dm e'^^^'^(™''l^)+'^"^\0)m (j 
{O}0{h) = Jdm e^[^Jv(A,/3)+fcm] ' • (29) 

The denominator is necessary because now the shifted effective potential is not nor- 
malised. 
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3. Description of the simulations 



The basic steps in a TMC simulation, which we shall discuss later, are the following. 

1. Select an appropriate sampling for M remembering that, essentially, M M + 
N/2. We shall discuss the choice of the m grid in Section 3.1. Naturally, wc must 
restrict ourselves to a finite simulation range, [mmin, Wmax], which introduces an 
(exponentially small) cutoff error. 

2. Run independent simulations for each m, measuring the tethered averages of h and 
the other relevant observables. 

3. The mean values {0)m,i3 are smooth functions of m, so they can be interpolated 
safely. Wc use cubic splines [31], but other methods may also work. Appendix A 
has some technical details about this and other points of our implementation. 

4. Integrate {h)rh,0 for the whole range of m. We use an average of the integral in 
both directions to reduce the systematic error: 



lN{m,f3) = X / dm' {h)^,^p - / dm' {h)r,,,^p . (30) 




This defines riAr(m, (3) up to an additive constant. Notice that this is not the same 
as forcing ^^{m^ [3) to be symmetric (which it cannot be, since m has a finite lower 
bound but extends to infinity). 

5. Normalise the pdf, 

/■'Tlmax 

c= / dm exp[Af/jv(TO,/3)]. (31) 

■'TOmin 

Then the effective potential is 

Qn {m,f3) = lN{m,f3)-j^ log c. (32) 

6. Compute the canonical averages with the interpolated tethered averages and equa- 
tion (28). The statistical errors are easily estimated with the jackknife method (see, 
e.g., [1]). The ith. block of our splines interpolates the ith jackknife blocks of each 
simulated tethered mean value. 

7. To obtain canonical results with a magnetic field h, simply substitute with i?^, 

f2%{m, p, h) = Qn {m, (3) + hm-j^ log c', (33) 

where c' is the new normalisation constant and use 

{O)0{h) = y dm e^<(™'^''^) (0)^,0. (34) 

Fig. 1 pictures this process. On the top panel we represent the energy density as a function 
of m, together with a horizontal line marking the canonical average at the critical point. 
On the bottom panel we plot the pdf p{rh). This function is reconstructed with great 
precision (we have to keep in mind that there is a factor N = Li^ vn the exponent). 
The rendered accuracy in pifn) allows us to obtain (u)/3^ = —1.419 05(5), correct to five 
significant figures, even though the range of variation of is of ~ 20%. The graphs 
were generated from the simulations described in Section 5.1. 

Once the values of m have been selected, the independent simulations arc carried out 
in a standard way. We use Metropolis dynamics to update the configuration. Let {(Jx) 
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Exact canonical value for L = 128 — 
Tethered values for L = 128 ^ 




Linear scale 
Log scale 



0.5 




Fig. 1. Computation of the canonical average from tethered averages at the critical temperature. The 
lower panel depicts the pdf of m, in logarithmic (left axis) and linear (right axis) scales, so that 
both the peaks and the tails can be seen. On the upper panel we plot the tethered mean values 
~("}m,/3c 3.S a function of m. The horizontal line is the exact value for L = 128 computed from [26], 
(u)^^ = —1.419 076... The continuous curves are our cubic spline interpolations (see Appendix). Our 
final result is (u)/3^ = —1.41905(5). Notice that the range of tethered values for (u)^,/3e several orders 
of magnitude greater than our statistical error. 

be the initial configuration and {cr^} the proposed new configuration (where one of the 
spins has been reversed). Then the probability of accepting the change is, from (23), 



0, 



{1 



e^S} 



if M' > M 
if M' < Af, 



where 



AS = -~(3{U' -U) + [M' -M) + 



N 
~2 



1 log 1 



M' - M 



(35) 



(36) 



M - M , 

An interesting point about our algorithm is that it is 'compulsory parallel'. The fact 
that we have to run simulations at several values of m adds a layer of trivial parallelisa- 
tion: we simply perform the simulations for each rh independently and only later combine 
their results to construct [2n (an operation that is essentially costless in computer time). 

A cautionary note: as our pseudorandom number generator we were originally using 
a 64 bit congruential generator reported in [32]. We simulated the lattice sequentially 
and extracted a random number per site, although in principle we would only need one 
when AS < in (36). This matches the conditions studied in [32,33], where significant 
deviations from the expected values were found using the same generator for a D — 3 Ising 
model. The error showed up only for large lattices, N > 128^. Perhaps unsurprisingly, we 
also obtained wrong results (farther than 3 standard deviations from the exact values) 
for a system with a similar number of spins {N = 1024^). This is our biggest system and 
we did not have any problems for all the smaller ones. Once this issue was identified, we 
added to the congruential generator a 64 bit version of the shift register method reported 
in [34] and redid all our computations. All the numerical results presented in this paper 
have been obtained using the sum of both generators (modulo 2^^). 
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3.1. How to select a good sampling of rh 

A good choice of the m we are going to simulate may reduce systematic and statistical 
errors significantly. Remember that while — l<m<l,min principle may extend to 
infinity. Actually, p{rh) has completely negligible values outside the range [—1/2,3/2], 
so we can restrict ourselves to that interval. If we look at Fig. 1 we see that p(rh) is a 
two peaked distribution, so values of rh inside the peaks contribute more that those in 
the middle or in the tails. These peaks get closer together and slightly narrower as wc 
increase L, so it may seem that the choice of m is quite delicate (specially considering we 
do not know p{rh) until we have run our simulation). ^ A different, but related, question 
also arises: is it better to compute; at more; points or more precisely with a coarser 

grid? We shall try to give some practical recipes below. 

The question is easy to analise for 7jv, see equation (30). Let us assume that we have 
obtained the mean value oi h &i K points in a grid. Our estimator [h]k is related to the 
actual value by 

[h]k^Wk + Tlk. k^l,...,K, (37) 

where rjk are the errors, expected to be Gaussian distributed, of similar size and statis- 
tically uncorrelated. Hence, our numerical estimate for 1^ will be given by a quadrature 
formula 

K K 

In ^'^gk{h)k + '^9kVk- (38) 

fc=i fc=i 

In this equation the first summand is subject to systematic error while the second one 
is the statistical error. Now, since the quadrature coefficients gk scale as 1/K, it is clear 
that the statistical error will scale as 1/VK. This suggests that doubling the number of 
points of the grid is equivalent to doubling the simulation time for each one. However, 
the analysis for canonical mean values is fairly more involved, since the errors in will 
be highly correlated for different m. Therefore, we perform a numerical experiment. 

Table 1 shows the values for the energy density at the critical point, — obtained 
in different series of runs. In the first column, we use evenly spaced points with 10^ 
Monte Carlo Sweeps (MCS) each. In the second column the points are also uniformly 
distributed, but now we perform 10'' MCS in each of them. The third and fourth columns 
are analogous, but with a greater density of points inside the peaks. We can reach several 
conclusions from this table: 

• If we use too low a number of points we will see a significant systematic error, no 
matter how precise they are. 

• Once the systematic error is under control, increasing the number of evenly distributed 
points has an effect of 1/^/ Appoints in the statistical error. This is best seen in the 
leftmost column. The effect is roughly the same if we increase the number of MCS in 
each point by the same factor (the errors of the first column are \/lO times greater 
than the corresponding ones of the second). 



* The discussion in this Section is relevant to the disordered phase and the critical region. In the broken 
symmetry phase, the peaks rapidly get very high and extremely narrow as we increase L and the criteria 
axe different (see Section 8). 
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Table 1 

Final value for —{u)fj^ as we change the number of points for the reconstruction of 17 jy a-nd their 
precision. MCS = Monte Carlo Sweeps for each point. The runs labelled 'uniform sampling' consist of 
.^^points values of m evenly distributed in the range [—0.4,1.4]. The runs labelled 'improved sampling' 
have 2/3iVpoints points evenly distributed in the same range, plus and additional Afpoi^ts/S inside the 
peaks, eflfectively doubling the density in the dominant regions. The last results of each column have a 
similar precision, but those computed with uniform sampling required twice the simulation time. 

Uriiforrii sampling Improved sampling 



JVpoints 10^ MCS 10'' MCS 10^ MCS lO'' MCS 



12 


1.43728(33) 


1.437 38(11) 






23 


1.419 25(22) 


1.419117(6) 






46 


1.419 21(13) 


1.419117(43) 


1.41908(11) 


1.419 107(38) 


91 


1.419 14(10) 


1.419 093(36) 


1.41913(8) 


1.419128(31) 


181 


1.419 14(7) 


1.419 095(28) 


1.419 034(5) 




451 


1.419 06(5) 




1.419 073(39) 




901 


1.419 065(33) 




1.419 077(26) 




1801 


1.419 062(24) 








Exact 




1.419 076 272 086. . . 





• If we add more points inside the peaks, the error may decrease faster than 1 / -^A^'pomts- 
The errors with iVpoints and uniform sampling are only slightly smaller than those with 
-^points/2 and improved sampling. 

We can summarise this analysis with the following prescription for the choice of m: 

1. Run a first simulation with enough uniformly sampled m so that the systematic 
error is small or unnoticeablc (i.e., so that the peaks of the distribution can be 
roughly reconstructed and the tails are reliably sampled). We have used ~ 50. This 
may seem a big number, but we must take into account that we have only looked at 
the energy in Table 1. Other quantities, such as high moments of the magnetisation 
or observables at a nonzero magnetic field, require that the tails of the distribution 
be reasonably well sampled. 

2. Add a similar number of points inside the peaks of the pdf to eliminate the sys- 
tematic error and reduce the statistical error. 

We have found that the second step is not always necessary. In fact, for lattices up to 
L = 256 we have limited ourselves to computing 51 evenly distributed m. For bigger 
lattices the peaks are steeper and we have added another 26 points inside them. 

3.2. Other practical recipes 

It is sometimes interesting to compute high moments of the magnetisation (for exam- 
ple, see Section 6). One obvious possibility is to simply measure individual values for 
m^(m; {(Ta;}) and then compute the tethered and canonical averages as usual. Btit our 
method provides an alternative way of calculating {m^)fj. Indeed, we have the whole pdf 
p{m) and we know that M = M + R. Now, the moments for R can be easily obtained 
analytically, so it suffices to compute (m^)/? to reconstruct {m^)p without any need for 
individual measurements of m^. For example. 
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{m^)p^{[n,-l/2f)p-^ . (39) 



These formulas are valid for the symmetric phase, where (m)/3 — 0. We have computed 
the moments of the magnetisation up to {rrfi) p both from individual measurements and 
with this procedure and the results are identical. This will not be at all surprising once 
we see Section 4, where it is shown that the correlation time for {m)m,i3 is < 1 (which 
means that the uncertainty in p(m) is going to determine the total error). 



4. Critical slowing down 



Before we examine our results, we must make sure that we have been able to thermalise 
our systems and that we have enough independent measurements to generate precise 
averages. We would also like to know whether our algorithm suffers from critical slowing 
down (CSD) and, if so, in what measure. 

To address these questions we have computed the autocorrelation functions and inte- 
grated autocorrelation times [5,1] for u, h, m and for each value of m (we will restrict 
ourselves to the critical temperature in this section). Another observable will be of inter- 
est. Remembering definitions (9) and (12), we can introduce the tethered mean value of 
the two point propagator at the smallest momentum: 




L = 


1024 1- 


L = 


512 ^ 


L = 


256 i- 


L = 


128 i- 


L = 


64 ^ 


L = 


32 1- 



Fig. 2. Integrated autocorrelation time for h and F as a function of rh, computed with the self consistent 
window method [5]. Notice the absence of critical slowing down for the former. In the lower panel, 
as the correlation time for F grows as L^, we did not use individual measurements of F after each 
Monte Carlo step for L > 512. Instead, we averaged over 50 such measurements for L = 512 and over 
4000 measurements for L = 1024. We then computed the autocorrelation functions for these blocked 
measurements (of course, we multiplied the integrated times thus obtained by the length of the blocks). 
This accounts for the smaller errors in Tint for these lattices. For L = 1024 the correlation time becomes 
unmeasurable (i.e., smaller than our blocks) for rh > 1.1 or m < —0.1. 
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0.01 



0.001 




Fig. 3. Normalised autocorrelation functions of u, F and h at the minimum (m = 0.51, black symbols) 
and one of the maxima (m = 1.14, grey symbols) of p{m) for L = 128 at the critical temperature. We 
appreciate how the curves for u (circles) and F (squares) become parallel after a small number of steps. 
We have only plotted pf^ (triangles) until the point where it first becomes compatible with zero, to avoid 
cluttering the graph. 



iki-x 



We define the autocorrelation function at time t for an observable by 

Co{t, m) 



Co{t,m) = {OsOs+t)A- (0) 



2 



po{t,m) = 



Co(0» 



From this function we can obtain several characteristic times [5 

t 



Tcxp,o("^) = lim sup ■ 



^og\po{t,m) 



' cxp 



sup Tcxp.O, 
O 



^ oo 

Tint,o{m) = X + X! m), 



exponential time of 0, 
exponential time of the system, 

integrated time of 0. 



(41) 

(42) 

(43) 
(44) 

(45) 



The first two measure the amount of time (i.e. number of Monte Carlo sweeps) that 
must pass before the system is thermalised. The third characterises the minimum time 
difference so that two measurements can be considered independent (i.e. uncorrelated). 

We compute Tint,o using the self consistent window method [5]. In Fig. 2 we show the 
integrated autocorrelation times for F and /? as a function of m. 

We see a clear difference: while the time scales for F grow as L^, with z w 2; the Tjnt 
for h do not exhibit any significant increase. That is, even though we are using a local 
algorithm, some of our observables do not experience any critical slowing down (and h is 
a particularly important observable, as we integrate it to obtain the effective potential). 
The integrated times for the energy are much smaller than those of F, but their behaviour 
is qualitatively similar (i.e., they grow as L^). Within the computed quantities, F is the 
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slowest mode for our algorithm. In some sense, the absence of CSD for all functions of m 
(including h) is not completely surprising, because the Metropolis update uses nonlocal 
information that involves the whole system through m. 

The exponential times are much more difficult to compute precisely, but we can easily 
give a rough estimate. Fig. 3 shows the normalised autocorrelation functions for the 
energy and F at two different values of rh (one at one of the maxima of p{rh), the other 
at the minimum). We see that F is a rather pure mode for our Metropolis dynamics, 
with a very clean exponential decay for small times, so that r^xp.F ~ Tint.F- The curve 
for the energy falls rapidly at first and then becomes parallel to pF{t,m). This is a clear 
indication that the integrated time for F can be considered as a good approximation to 
the exponential (i.e. thcrmalisation) time of the system. 

The reader may find this section in conflict with the exact result by Hohenberg and 
Halperin [35], who showed that for conserved order parameter dynamics (model B) in 
Ising models the critical exponent is 76 = 4 — r/. Nevertheless, even our slowest mode has a 
much smaller z. The way out of this paradox is in the nonlocal nature of our conservation 
law. Model B needs magnetisation diffusion across the boundary if the magnetisation 
is to change in the enclosed region. When the locality constraint is violated, the new 
dynamical exponent is expected to be Znoniocai = z — 2 [36]. Clearly, this is the case for 
our dynamics, where a change of rh inside a given region does not occur through diffusion 
across its boundary. 

5. Results at /9c ) h = Q 

We summarise here our results at the critical point, which correspond to the bulk of 
our simulations. This section is divided in three parts. First we provide the parameters 
of our simulations. Second, we report our results for the mean values that would be 
considered in a traditional canonical computation. Finally, in subsection 5.3 we perform 
a rather unconventional, but very accurate, computation of the critical exponents ratio 
(i/v taking advantage of the peculiarities of the tethered formalism. 

5.1. Parameters of our simulations 

We have simulated lattice sizes L = 16, 32, . . . , 1024. For systems with L < 256 we have 
used 51 uniformly distributed values of m in [—0.5, 1.5] (except for X = 16, where we had 
to extend the range to [—1.61, 1.61] to avoid cutoff errors). As we shall see in Table 2, this 
choice (considered suboptimal in Section 3) does not introduce any discernible systematic 
error. We have been very conservative and chosen a wider range for m than what we used 
in Section 3.1, slightly enlarging the computations to avoid any chance of a cutoff error. 
For L > 512 we ran a first simulation with 51 equally spaced m in [—0.3, 1.3]. In these 
cases the peaks were considerable narrower and closer together than for the smaller 
lattices. To avoid any discretisation errors, we have followed the practical recipe given in 
Section 3: after obtaining an approximation to the shape of J7jv, we have added another 
26 points, doubling the density inside the peaks. 

In all cases we have performed 10^ Metropolis sweeps for each value of m. Following the 
standard prescription [5] , we have discarded a fifth of the measurements for thcrmalisation 
(although the correlation times are much smaller, see Section 4) and formed 100 jackknife 
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blocks [1] with the rest to obtain a rehable estimate of the statistical error. According 
to Section 4 and, particularly, Fig. 2, this means that for L = 1024 the length of the 
blocks is only about 5 times the largest found Texp (that of to = 0.5). There has been no 
need to increment the numerical effort in this point as this exponential time corresponds 
to the minimum of p(to) and wc have seen that the canonical values are dominated by 
the neighbourhood of the peaks. There the exponential time is an order of magnitude 
smaller, so our error estimates are sound (as can be seen by comparing our results with 
the exact values). In any case, we have recomputed the errors of Table 2 with 50 and 
200 jackknife blocks and checked that they remain completely stable. Our canonical 
simulations consisted of 10^ Swendsen-Wang updates. 

5.2. Canonical averages at the critical temperature and zero field 

Wc present in Table 2 our values for the canonical averages of several physical quan- 
tities. We have compared with the exact results for and C in finite lattices computed 

Table 2 

Results at the critical temperature and comparison with a cluster algorithm. (T): Tethered Monte Carlo, 
(C): Swendsen-Wang, (E): Exact results at finite L from [26]. 



L 


-<">/3c 


X2/L^ 






C 


B 


a/3€i/io4 


16 (T) 


1.453 08(4) 


0.545 43(6) 


0.9116(2) 


0.246 13(13) 


7.718 6(14) 


1.165 62(7) 


0.036 547(9) 


16 (C) 


1.452 9(2) 


0.5451(3) 


0.910 4(9) 




7.718(10) 


1.165 9(3) 


0.036 50(6) 


16 (E) 


1.453 065... 








7.717134... 






32 (T) 


1.433 69(4) 


0.459 00(10) 0.9072(4) 


0.242 2(3) 


9.509(3) 


1.16723(14) 0.144 07(7) 


32 (C) 


1.433 67(12) 


0.459 1(2) 


0.9078(9) 




9.493(13) 


1.1671(3) 


0.144 1(3) 


32 (E) 


1.433 659. .. 








9.509 379... 






64 (T) 


1.423 97(4) 


0.38619(18) 


0.906 5(9) 


0.240 0(5) 


11.285(6) 


1.1675(3) 


0.573 8(6) 


64 (C) 


1.423 90(6) 


0.386 0(2) 


0.905 6(10) 




11.293(17) 


1.167 7(4) 


0.5731(11) 


64 (E) 


1.423 938... 








11.288 138. . 






128 (T) 


1.419 05(5) 


0.3244(3) 


0.904 0(18) 0.240 8(11) 


13.063(10) 


1.168 4(7) 


2.289(5) 


128 (C) 


1.419 06(4) 


0.324 59(17) 0.904 8(10) 




13.06(2) 


1.1677(4) 


2.287(4) 


128 (E) 


1.419 076. .. 








13.060 079. . 






256 (T) 


1.416 63(5) 


0.272 8(6) 


0.904(4) 


0.240(2) 


14.83(2) 


1.168 7(14) 


9.16(4) 


256 (C) 


1.416 64(2) 


0.272 86(14) 0.904 2(9) 




14.83(2) 


1.168 2(4) 


9.14(2) 


256 (E) 


1.416 645... 








14.828 595. . 






512 (T) 


1.415 42(4) 


0.229 3(7) 


0.903(6) 


0.240(4) 


16.57(3) 


1.168(2) 


36.38(19) 


512 (C) 


1.415 444(11) 


0.229 68(13) 0.905 9(10) 




16.60(2) 


1.1676(4) 


36.64(10) 


512 (E) 


1.415 429... 








16.595 404. . 






1024 (T) 


1.414 89(4) 


0.194 9(15) 


0.919(15) 


0.240(10) 


18.28(8) 


1.163(6) 


148.8(19) 


1024 (C) 


1.414 826(6) 


0.193 07(12) 0.904 6(11) 




18.35(3) 


1.1681(4) 


146.1(4) 


1024 (E) 


1.414821. .. 








18.361348. . 
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with the expressions given by Ferdinand and Fisher [26]. We have also run simulations 
with a Swendsen-Wang cluster algorithm (we use our own implementation, based on the 
one distributed with [1], but our results are compatible with those of [37]). 

From Table 2 we can confirm that TMC is capable of producing very accurate results. 
The relative errors for X2 and B scale as L. This can be accounted for by noticing 
that both are completely determined by p{rn), recall Section 3.2. Indeed, while h is self 
averaging and virtually free of critical slowing down (meaning that for a fixed simulation 
length its error scales as l/\/]V), we are multiplying {2^(171,(3) by a factor of A^, yielding 
an overall ^/N scaling for the errors. 

As we said in the Introduction, TMC is not meant to be a competitor to cluster 
algorithms for the Ising model without magnetic field. For example, the CPU time to 
compute each of the 77 simulations at fixed rh for L — 1024 is about a third of what we 
needed for the whole Swendsen-Wang simulation, which is also more precise. 



5.3. Finite Size Scaling and the peaks of p\{m) 



Let us obtain an accurate estimate of the critical exponent ratio P/v (known to be 
1 /8) in a way that would not be competitive for a canonical computation. 

Our starting point is the Finite Size Scaling formula for an arbitrary observable O 
(see, e.g., [1]), as we get closer to the critical point {(3 = (3^, h = 0): 

{0)t{h)=L-=^°/''\fo{L^^''t,Ly-h) + ..], i=^V^- (46) 

Here the dots represent possible corrections to scaling. We shall work at = 0, so we 
only have to consider the first argument of the scaling function fo- 

Elaborating equation (46) and recalling that /3 is the critical exponent for the mag- 
netisation, it can be shown that (sec, e.g., [1]) 

p,{m,p,;L) = L^/''f{L''/''m), (47) 

where pi{m, (3c', L) is some smoothed version of pi(m, /3c; i), equation (16). Recalling 
that p{m, pc',L), equation (20), is just one of such smoothings, we can substitute this 
expression by 

p{m, f3c; L) = L^/"" f {L^''' [m - 1/2)). (48) 
In particular, the pdf has two maxima at nn^^^^^ + i, whose scaling behaviour is expected 
to be m± ,k <x 

Recall that we measure directly /?, which is the derivative of the logarithm of p(to) and 
thus is exactly zero at these peaks. Therefore, for each of the maxima we just have to 
identify the two consecutive points of the grid such that {h)rhi,0 > and {h)rhi+x,l3 < 
and find the root of the cubic spline joining them (we have also used jackknife blocks to 
estimate the errors). 

The position of these peaks is a canonical observable, but one that is more easily 
obtained through the tethered ensemble (in a canonical simulation we would have had to 
construct pi{m) and locate its maximum directly, a harder problem than finding a zero). 

The results of this analysis for the simulations described in Section 5.1 are collected 
in Table 3. Notice that even with our relatively coarse rh sampling we have been able to 
determine the position of the peaks with a precision ranging from 0.013 % to 0.46 %. If 
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we had wanted to give a very precise value for these points, we could have done so with a 
rather small investment in CPU time, as we only need precise simulations at two values 

of rh, conveniently chosen. 

By fitting 'rn^^^y.{L) to a power law we can obtain an estimate of (i/u. Table 4 gathers 
the results of fitting the data of Table 3 to a power law, for lattice sizes L > i,nin- We 
see that for Ly^in — 32 the power law does not represent the curve adequately. To give 
an error estimate that represents both systematical and statistical sources we follow the 
criterion explained in [38]: when two consecutive values are first compatible, wc give the 
most precise as central value, but with the bigger error. For the negative magnetisation 
peak we find p/u = 0.123 9(11) and for the positive peak f3/u = 0.124 5(10). 

Following [37] , wc can even go further and try to characterise the corrections to scaling. 
We assume that in the Ising model in Z? = 2 the dominant corrections to scaling are 
analytical 

m%^^ = L-^/''[A^+B^L-% A = 7/4. (49) 

We have fitted our points for all lattices to this expression, fixing the exponents to their 
exact values and varying A"^ and B^. We have obtained x^/d-o-f- = 0.9858/4 for the 
negative peak and x^/d.o.f. = 2.825/4 for the positive one. 

Table 3 

Position of the peaJcs of the probability density function of the magnetisation for several lattice sizes. 
As we are at the critical point, the values in the table extrapolate to zero when L — > oo (compare with 
Section 8, where we study the ferromagnetic region). 



32 0.764 01(10) 0.764 31(11) 

64 0.702 86(18) 0.703 0(2) 

128 0.645 3(3) 0.645 1(4) 

256 0.5921(7) 0.5910(7) 

512 0.5419(12) 0.542 7(9) 

1024 0.499(2) 0.500(2) 



Table 4 

Fits of mpg^^(L) to a power law, "ipg^^j. = ^, including all lattice sizes L > Lmin- We give the 
computed exponent and the chi square per degree of freedom. As discussed in the text, for small lattices 
we detect corrections to scaling. Our fits converge to the exact result, /3/i/ = 0.125. 



peak 



'min 


/3/I/ 


xVd.o.f. 


13/u 


xVd.o.f. 


32 


0.1217(3) 


23.44/4 


0.122 4(3) 


27.85/4 


64 


0.123 9(5) 


2.027/3 


0.124 5(5) 


2.087/3 


128 


0.125 0(11) 


0.7569/2 


0.124 6(10) 


2.053/2 


256 


0.126(4) 


0.6456/1 


0.122 0(23) 


0.3248/1 
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6. The scaling paramagnetic region: the renormalised coupling constants 



We shall show here that TMC works as well in the paramagnetic phase in the scaling 
region. For definiteness, we shall study the renormalised coupling constants (see [2,3]). 
These constants arc notoriously difficult to compute using Monte Carlo methods, so they 
provide a good challenge for our formalism. 

Let us consider the D = 2 Ising model in its thermodynamical limit. Then we can 
define the Gibbs free energy by 

Git, h) = hm 1 log[Z(i, h)l t= (50) 

and its Legendre transform, the Hclmholtz free energy (mj = {m)t), 

F{t,mt) =max[mth-G{t,h)]. (51) 

From a field theoretical point of view, the latter can be expanded as a series in the 
renormalised coupling constants g2n (here we shall use the definitions of [39]) 

Fit, mt) - Fit, 0) = 1 f + f: ^<P'-] (52) 

(53) 

In these equations. 

We would like to express the g2n in a form suitable for our lattice simulations. To do this 
we remember definition (15), according to which Git, h) can be immediately represented 
as a Taylor expansion with coefficients X2n, 

00 

Git,h)-Git,Q) = J27^h'''- (55) 




Combining all these equations we arrive at the following explicit formulas: 

X4it) 

t-'?+X2(i)W 



54 = - lim ^^1^, 52n = r-2„(ff4)"-\ (56) 



with 



,, = 10- lim ^^^1^, (57) 



t- 



rs = 280 + lim 



X8(t)X2(t)' _ggX6(0X2(0 



X4{tr xiity 



(58) 



There has been a great deal of interest about these coupling constants and many precise 
computations have been performed, both with field theoretical and numerical methods. 
Balog et al. [40] arrive at §4 = 14.697 5(1) with a form factor approach while Caselle 
et al. [39] give 54 = 14.697323(20), re = 3.678 66(3)(2) and rg = 26.041(8)(3) using 
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transfer matrix techniques. There are fewer Monte Carlo determinations of these quan- 
tities [40,41]. 

In this Section we shall reproduce the Monte Carlo computations of [41] with the 
TMC method to give our own estimate of the first coupling constants. A Monte Carlo 

determination of the g2n can obviously not be directly computed for an infinite lattice. 
Instead, we have to perform a double limit limt_^o+ liniL-^oo, that is, run simulations very 
close to, but above, the critical temperature at increasing values of L until the results 
are stable (within our errors). Ref. [41] concludes that L = 100 is a big enough lattice. 
In order to compare our results to those of that reference, we have worked at the same 
temperature, /3 = 0.42. 

We have used 191 uniformly spaced m in [-0.45, 1.45], with lO'^ Monte Carlo sweeps 
each. Table 5 summarises our rcsiilts, which arc compatible with those of [41] but more 
precise. These calculations involve working with high moments of the magnetisation (we 
have used the indirect method explained in Section 3.2, combining the cumiilants of r 
and m). The table also compares the tethered results with a Swendscn-Wang simulation 
(SW). Notice that in the latter the errors increase very quickly when we consider high 
powers of m. For instance, for L = 100 the error in the susceptibility is almost the 
same for the TMC and the SW methods, while for gs the error of the latter is almost 
ten times higher. Our whole TMC simulation for L = 100 required about 100 hours 
of computer time, while the Swendscn-Wang one was completed in about 10. If we use 
'improved' or cluster estimators (SWC) [42], however, the advantage of the tethered 
algorithm disappears. 

Table 5 

The renormalised coupling constants. TMC: Tethered Monte Carlo; SW: Swendsen-Wang; SWC: 
Swendsen-Wang with cluster estimators. For L = 100 we also give the results of Ref. [41], computed 
with a single cluster Monte Carlo method. 







L = 50 






L = 


100 






TMC 


SW 


SWC 


TMC 


SW 


SWC 


Ref. [41] 




1.228 238(14) 


1.228 31(4) 


1.226 067(7) 


1.226 076(8) 




X2 


196.85(9) 


196.96(15) 


196.91(5) 


203.78(11) 


204.07(10) 


203.92(2) 


204.4(3) 


«1 


11.749(5) 


11.756(9) 


11.753(3) 


11.888(11) 


11.907(10) 


11.893 2(10) 


11.90(2) 


re 


4.4462(9) 


4.449(4) 


4.446 9(10) 


3.70(6) 


3.73(8) 


3.731(6) 




ra 


39.76(3) 


39.83(13) 


39.77(3) 


26.2(6) 


24(3) 


26.47(18) 




94 


12.817(6) 


12.80(3) 


12.812(7) 


14.66(5) 


14.69(9) 


14.673(8) 


14.60(16) 


96 


730.4(5) 


729(3) 


729.9(9) 


794(9) 


806(24) 


803.3(16) 


8.5(4) • 10^ 


98/10^ 


8.370(8) 


8.36(6) 


8.363(16) 


8.25(13) 


7.5(11) 


8.34(7) 


8.8(10) 



7. Results at j3c, h ^ Q 

An interesting feature of TMC is that it allows us to obtain accurate data in the 
presence of a magnetic field. In order to do this, we reanalise the data from the simulations 
described in Section 5.1. Let us stress the fact that we do not have to run any new 
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simulations at all, we simply use the tethered values computed at zero magnetic field 

and modify the effective potential as in equations (33) and (34). 

We will typically be interested not in computing the equivalent of Table 2 for a par- 
ticular value of h, but in drawing a curve {0)/}{h). When doing this, we can improve our 
statistical errors somewhat if wc take into account that the observables can be either odd 
or even in h. This means that wc can (anti)symmetrise the curves: 



iOrrih) = iOM^)-^iOM-h) , (59) 

Now, if the data were completely uncorrelated, so that {0)p{±h) were independent, this 
operation would imply a l/-\/2 reduction in the statistical error. Actually, this is not the 
case: the individual values for ±h are very strongly correlated and the symmetrisation 
reduces only very slightly the error for even quantities. For odd observables, however, 
the fact that we are computing a difference rather than a sum yields a very significant 
decrease in the statistical error (around a factor of 10, specially for small values of h). 
Wc shall use equations (59) and (60), but wc shall drop the explicit 'odd' or 'even' 
superscripts. 

The first observable wc shall consider is the correlation length £,2{h), defined in equa- 
tion (13). This is even in h, so we can use the symmetrised version of equation (60). 
In this case, unlike previous sections, no data for ^2(^) arc readily available in a finite 
lattice and doing our own cluster simulations would be harder (as we pointed out in 
the Introduction, the currently most popular methods in the h ^ regime are transfer 
matrix techniques). Nevertheless, we can provide a very good check of consistency using 
the FSS equation (46). Indeed, given that our simulations have been carried out at the 
critical temperature, we can consider the critical behaviour as /i — > 0. In our case, we are 
working at precisely /3 = /3c, so the first variable in the scaling function disappears and 
we have {xq = for the correlation length) 

^2/L^f^,{Ly^h). (61) 

The critical exponent yh is 15/8 for the two dimensional Ising model and the function 
is expected to be very smooth. As we can see in Fig. 4, equation (61) is perfectly 
valid in our case, once we discard the data for L < 32. At a first glance, it may seem that 
we have even overestimated our errors for L = 1024, but remember that the points are 
very strongly correlated. The universal scaling curve is well represented by a sixth order 
even polynomial, with a value of diagonal x^/d.o.f. of 3.978/36.^ 

Thus, the universal scaling function /^2(a;) for x < 1.5 is extremely well represented 
by 6 

f^2ix) = 0'0 + a2x'^ + a4x'^ + aQX^, (62) 



^' The points of the graph are very strongly correlated, so their covariance matrix is singular and the 
parameter of the fit cannot be computed. Instead, we restrict ourselves to the diagonal elements 
of the covariance matrix and minimise the resulting 'diagonal' Xd- Usually, the degrees of freedom are 
computed as number of points minus the number of parameters. In our case this does not correspond to 
the actual degrees of freedom of the fit, but we have maintained the usual notation. 
^ To estimate the errors we have computed a fit for each jackknife block. 
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Fig. 4. Finite Size Scaling with (2{h). We represent the scaling function by a fit to 
f^{x) = ao + a2X^ + 042;* + afjX^. The value of the diagonal chi square per degree of freedom is 
X^/d.o.f = 3.978/36, which confirms that the points fall on a smooth curve within their errors. 

with 

ao = 0.239 9(2), 03 = -0.063 9(4), 04 = 0.023 5(7), ag = -0.004 5(4). (63) 

We now consider the magnetisation, m^^(/i) ~ {m)p^{h). Contrary to the correlation 
length, the magnetisation is an odd function of h. This means that by antisymmetrising 
our results for ±ft, as in equation (59) we can greatly reduce the error. In this case we 
again lack a convenient way of generating the same curve with a different method, but 
we still can check the validity of our results. To do this, we differentiate equation (55) 
and notice that the coefficients of the Taylor expansion of {m)p^{h) are none other than 
the 2n point correlation functions at zero momentum: 

m0^{h)^{m)p^{h)^——^=X2h+^h + —h + —h +.... (64) 

This expansion provides a new way of computing X2n- we just generate a reasonable 
number of points of the m(/i) curve, which we parameterise with a truncated version 
of equation (64). The choice of values for h is somewhat delicate: if we use very small 
magnetic fields we will only be able to appreciate the first few coefficients but if we go 
too far in h wc would need to have sampled the tails of the pdf of to very precisely. We 
have found that magnetic fields up to ^ (X2)^^ provide a good compromise. 

A second difficulty is that posed by the correlations among the points of the curve. 
These are so strong that the covariance matrix turns out singular, which bars us from 
obtaining a fit and its errors by minimising x^- Instead, we have computed an odd 
interpolating polynomial with a finite difference formula (this gives us as many X2n a-s 
points) and we have tried to control the correlation by estimating the errors with the 
jackknife method. This takes care of the statistical error. To control the systematic error, 
we have both reduced the range in h and varied the number of points to check whether 
the parameters were stable. We have found that if we use this method with n points 
the last one or two coefficients are usually unreliable (i.e., their value changes beyond 
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Nonlinear susccjjtibilitics from direct iiicasureriiciits at /?, = 


(M) and from a, finite 


difference formula 


foi' /// ;^ ih) 


as a fund ion of llie nia,i^nel 


ic field (F). 






L 






N-^XO 




16 (M) 


0.545 43(6) 


-0.545 72(13) 


2.265 7(8) 


-20.059(10) 


16 (F) 


0.545 43(6) 


-0.545 7(2) 


2.262 8(19) 


-19.70(14) 


32 (M) 


0.459 00(10) 


-0.386 1(2) 


1.348 5(10) 


-10.042(10) 


32 (F) 


0.459 00(10) 


-0.386 2(3) 


1.3484(18) 


-10.00(2) 


64 (M) 


0.386 19(18) 


-0.273 3(3) 


0.8031(13) 


-5.032(11) 


64 (F) 


0.386 19(18) 


-0.273 8(5) 


0.805(2) 


-5.02(2) 


128 (M) 


0.3244(3) 


-0.192 8(4) 


0.475 8(17) 


-2.504(12) 


128 (F) 


0.324 4(3) 


-0.194 3(7) 


0.481(3) 


-2.52(2) 


256 (M) 


0.272 8(6) 


-0.136 2(7) 


0.283(2) 


-1.250(12) 


256 (F) 


0.272 5(6) 


-0.135 8(15) 


0.280(6) 


-1.20(4) 


512 (M) 


0.229 3(7) 


-0.0964(7) 


0.168(2) 


-0.625(9) 


512 (F) 


0.229 3(7) 


-0.096 0(12) 


0.166(4) 


-0.60(2) 


1024 (M) 


0.194 9(15) 


-0.069 8(13) 


0.104(3) 


-0.328(12) 


1024 (F) 


0.194(3) 


-0.063(6) 


0.08(2) 


-0.2(3) 



the error bars if we add another point). This only means that if we want to obtain n 

physically meaningful parameters, wc should compute at least n + 2 points. Wc want to 
compute the nonlinear susceptibilities up to xs, so to be safe we have used 7 points for 
each lattice size, equally spaced at intervals of (10x2)~^, where X2 is the susceptibility 
computed from the simulation at h = 0. 

Table 6 summarises our results for the X2m computed from the measurements at /i = 
and from the finite difference formula a.t h ^ 0. We see that both series of values are 
compatible, but that the former arc more precise. 

We can also do a FSS analysis with the magnetisation. Notice that at small applied 
fields m0{h) ~ x^h. If we increase the field, however, we can appreciate a deviation from 
the linear behaviour (sec the left graph of Fig. 5) . If wc want to collapse all the curves 
on one we must take the critical exponent of the magnetisation, — /3, into account: 

mp^{h)c^L~^/''f^iLy'^h). (65) 

The function fm is depicted on the right panel of Fig. 5. We fit it to an odd polynomial 
fm{x) = aix + a^x^ + a^x^ + arx"^ for lattices L > 64, as we did for the correlation 
length. The value of the diagonal Xd/d-O.f. for this fit is 41.85/31. The last point of the 
curve, which corresponds to the highest magnetic field for L = 1024, seems to show a 
small deviation from the curve. The reason may be that we are taking the range of the 
scaling variable too far and that scaling corrections start to act. Remember that wc had 
spaced the values of h in units of (10x2)~^, so that the representation on the left scaled 
properly, which is not be best choice for the FSS analysis. 

It is interesting to examine the behaviour of the magnetisation with small but L 
independent magnetic fields (Fig. 6). We observe two well differentiated scales: an FSS 
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Fig. 5. Magnetisation at nonzero magnetic field for different lattice sizes. We provide two representations. 
The one on the left is plotted in terms of the applied field times the susceptibility at h = 0. This way 
we can appreciate the linear and nonlinear regimes. The graph on the right is the FSS curve for m^(/i) 
(a fit to a seventh order odd polynomial, with diagonal Xj/d.o.f = 41.85/31). 



0.8 - 




Fig. 6. Magnetisation at the critical point for several lattices and an L independent range of magnetic 
field values (notice the logarithmic scale of the horizontal axis). 

regime, where the slope of the curve is very large, and a thermodynamic limit, where the 
curves merge. The range of h is limited by the largest obtained {h)m,i3 in the simulated rfi 
grid (this is the reason for the rapid growth of the error bars for the final points, specially 
noticeable in the largest lattices). Due to the small value of P/iy, the displacement of the 
curves for small h is almost linear in logL, see equation (65). 

We can repeat the process described here for other observables, such as the energy 
or specific heat. However, our computation can be inefficient if we do not prepare it 
carefully. Remember that for h — the canonical average was dominated by two peaks, 
whose position was determined by the value of rh such that {h)^ = 0. Now the value 
with a nonzero magnetic field will be even more dependent on a saddle point, determined 

by 

{h),h = h. (66) 
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In practice, this means that the way to obtain a precise result for, say, the energy at a 
given field h is to first estimate the value of m such that the previous equation is satisfied 
and then run a long simulation there. 

8. The broken symmetry phase 

Let us now consider the (3 > f3c regime. In this situation, the infinite system shows 
a nonzero expectation value for the order parameter, Ta^^fj^ = {m)f}yf3^ ^ 0, even in 
the absence of an external magnetic field. This may seem incompatible with the parti- 
tion function (1), where the configurations {ctcc} and occur with equal probability. 
The well known solution for this apparent paradox is spontaneous symmetry breaking [2] , 
whose mathematical formulation involves considering a small magnetic field (which es- 
tablishes a preferred direction) and taking the double limit 

(m)/3.oo = lim lim {m)i3x- (67) 

The order of the two limits is crucial: were we to reverse it the magnetisation would 
always vanish. We see then that the symmetry of our model complicates the definition 
of a broken symmetry phase for finite lattices in the canonical ensemble. The traditional 
workaround consists in considering not the magnetisation m, but its absolute value \m\. 

The tethered ensemble provides a cleaner concept of broken symmetry phase. Consider 
the pdf of m, as in Fig. 1. In the ferromagnetic phase the corresponding graph will again 
have two peaks, but now these will be much narrower and higher, approaching two 
Dirac deltas in the thermodynamic limit. Suppose we want to perform the double limit 
of equation (67). This would involve introducing a small magnetic field which would 
shift the origin of h (29). The neighbourhood of one of the peaks would then become 
exponentially suppressed and eventually disappear in the thermodynamic limit. Thus, 
we can mimic the effect of equation (67) by considering only one of the two peaks from 
the beginning. This would not work below Pc, as there the peaks extrapolate to m = 
(recall Section 5.3). In a more complex model, the L evolution of the right peak should 
be checked, in order not to mistake a paramagnetic phase for a ferromagnetic one. This 
criterion suggests that TMC could be a powerful method for determining the order of 
magnetic phase transitions. 

This procedure has the considerable advantage that it works for any lattice size. In 
this section we have chosen the peak of positive magnetisation. We will illustrate it by 
considering the thermodynamic limit in the canonical ensemble in Section 8.1 and by 
studying the equivalence between the tethered and canonical ensembles in Section 8.2. 

8.1. The thermodynamic limit 

We have run simulations for lattice sizes L = 128,256,512, 1024 at /? = 0.4473 > /3c- 
This temperature was chosen because we estimated that it would roughly mirror the 
value of the correlation length ^ for our paramagnetic simulations. 

Following the previous discussion, we have worked in the rh > 0.5 (positive magneti- 
sation) region, where there is only one peak. An appropriate sampling of rh is even more 

In the ferromagnetic regime is not a good definition and we always use ^2, see equation (13). 
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Fig. 7. Computation at /9 = 0.4473 > /3c (ferromagnetic phase). The peak of the pdf of m gets narrower 
and closer to my + 1/2 as L increases (the vertical line, where my is Yang's magnetisation for the infinite 
system at /3 = 0.4473 [43]). Compare the scale of the x axis and the height of the peaks with those of 
Fig. 1. 



important in this phase (but easier to optimise) than in the situation described in detail 
in Section 3.1. The reason is that the peak is now so narrow that a choice of m spaced 
as in the aforementioned section would not only be completely wasteful, but may also 
completely fail to sample the peak. 

In the case of the Ising model, we know Yang's exact solution mY(/3) for the mag- 
netisation of the infinite system [43]. The positive peak for p{m) will then be very close 
to mY(/3) + ^ and get closer as we increase L. With this information in hand, we can 
adequately reconstruct the effective potential by running simulations in a small neigh- 
bourhood of myiP) + 5. For a different model, where we would lack the knowledge of the 
peak's position in the thermodynamic limit, we can just run simulations with a very fine 
grid for some small and essentially costless lattice size and infer from them an efficient 
distribution of points for the larger systems. 

We have represented p{m) for all the simulated lattices in Fig. 7, which includes the 
whole range of m for L > 512 (for the smaller lattices we have used a somewhat larger 
interval). It is interesting to compare the scale on the x axis with that of Fig. 1. As will 
be discussed in detail in Section 8.2, the peak approaches mY(/3) + ^ (the vertical line) as 
L increases. Table 7 compares the values of the energy and specific heat obtained in our 
simulations with the exact values given in [26]. Notice how very small simulated ranges 
of m {Am — TOmax ~ ?Timin) yield very accurate results. In fact, we can see that for the 
L = 1024 lattice we obtain a more precise determination for the energy with 27 points 
than what we obtained at the critical temperature with 77 (we still perform 10^ Monte 
Carlo sweeps in each point). This result is even more impressive if we consider that some 
of these 27 points, being deeply inside the tails of the distribution, do not have any effect 
whatsoever in the average with our error (of course, we do not know this until we have 
run the simulation and seen the actual width of the peak). 

From Table 7 we can conclude that the thermodynamic limit has already been reached 
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for L = 256, at least to the level indicated by our errors. Our whole computation for 

L = 512 required about 270 hours of computer time. For comparison, a 30 hour long 
Swendsen-Wang computation of the correlation length for L = 512 gives ^2 = 11.8(2). 
We see that the ratio of computation time for both methods has changed significantly 
from the critical point, where the advantage of the cluster algorithm was much greater. 

Let us now consider the curve for m^(/i) = {m)p(h) in the thermodynamic limit. We 
compute it with the same method employed in Section 7, but now we cannot apply 
the antisymmctrisation (59). The result is plotted in Fig. 8, where we plot L = 512 for 
magnetic fields in the range h E [0, 10^'^]. We also plot L = 256 to check that both curves 
coincide and we have reached the thermodynamic limit within our errors. We appreciate 
in this figTire just how efficient the antisymmctrisation procedure was in Fig. 5, which 
had much smaller errors. 
Table 7 

Canonical averages for several physical quantities of an Ising lattice at /3 = 0.4473 computed with the 
tethered method (T). The grid of m values is uniform in the narrow simulated band. Also included are the 

exact results for finite lattices from [26] and the exact results in the thermodynamic limit from [12,43]. 
We appreciate that by simulating only a very small range Am of values for m we can obtain very precise 
values. Within our error, we have already reached the thermodynamic limit for L = 512. 



L 


^points 


Am 


-{u}0 


C 


€2 






X2 


128 (T) 


90 


0.9 


1.490 397(18) 


8.874(4) 


10.394(17) 


0.519 87(8) 


0.719 34(6) 


39.65(8) 


128 (E) 






1.490 409 763. . . 


8.877 363. . . 










256 (T) 


79 


0.39 


1.490 407(11) 


8.869(5) 


11.26(4) 


0.518 16(5) 


0.719 41(4) 


39.36(7) 


256 (E) 






1.490415 672... 


8.874 075... 










512 (T) 


27 


0.13 


1.490419(5) 


8.877(5) 


11.5(3) 


0.51777(4) 


0.719 45(3) 


39.37(9) 


512 (E) 






1.490415 689... 


8.874 046... 










1024 (T) 


27 


0.13 


1.490 416(4) 


8.868(7) 


11.4(18) 


0.51764(4) 


0.719 45(2) 


39.48(11) 


1024 (E) 






1.490 415 689. . . 


8.874 046. .. 










00 (E) 






1.490 415 689. . . 


8.874 046. .. 






0.719 436... 





8.2. Ensemble equivalence in the ferromagnetic phase 

Once we wander away from the critical point, the main objective is finding the value 

of physical quantities in the thermodynamic limit. The ensemble equivalence property 
suggests a way to reach this limit without constructing the whole canonical p{iri), but 
by concentrating instead on its maximum. Prom the computational point of view, this 
supposes a dramatic reduction in the needed effort for a TIVIC simulation. 
Ensemble equivalence can be expressed in mathematical terms by 

hm {0)f, = lim (0)(rf,),,^. (68) 

This equation can be understood as a more formal way of summarising the behaviour of 
Fig. 7. Indeed, we saw in the previous section that we could reconstruct the canonical 
averages considering only a very narrow range of m; in the thermodynamic limit a single 
point would be sufficient. 
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Fig. 8. Magnetisation in the presence of a magnetic field in the ferromagnetic regime for /3 = 0.4473. 
The horizontal line is the exactly known thermodynamic limit for h = 0. 

For the Ising model we can situate this point exactly because we know Yang's sponta- 
neous magnetisation 

my(/3) = hm {m)p = lim {m)p + ^ = [l - (sinh 2/3)"^] + i (69) 

We could then run simulations for several lattice sizes at precisely rhy and study the 
evolution of {0)my,i3 as we increase L. This is not the most practical approach, as for a 
model other than the D = 2 Ising lattice we would not know the position of the peak 
beforehand. Instead, we will follow a more general analysis that would work in more 
complex situations. 

Let us consider the canonical average of some quantity and recall that we are using 
periodic boundary conditions, so the approach to the thermodynamic limit is exponential 



(0) 



13 



1/2 



dm p(m, L){0),n = Of + AqQ 



(70) 



where Aq is a constant amplitude. We have just considered the positive magnetisation 
peak. 

Since the integral will be dominated by a saddle point at rh^^^^, with Wpoak ttiyi 
we can approximate the pdf by a Gaussian 



p(to,/3; L) 



! N 



■ exp 



N{m - rripcak)" 



2X2 



(71) 



Therefore, we expect the tethered average of at this saddle point to approach the 
canonical average (70), with a correction of order N^'^, 



Of = {0)p - Aoc 



To ease the notation we shall use the definition 



(72) 
(73) 
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Table 8 

Tethered mean values of several parameters at the peak of the probability density function for (3 = 
0.4473 and j5 — 0.6, together with the value of = {h)mY (this observable is zero at the peak and helps 
characterise how close we are to it). The exact value for an infinite lattice, which coincides with the 
canonical average, is also included for comparison. 







/3 = 0.4473 






/3 = 0.6 




L 


■ 10^ 


'^peak 


'"■peak,/3 




'^^eak 


■"peak,/3 


16 


2984(3) 


1.343 84(4) 


1.577 07(8) 


-290.6(14) 


1.471943(7) 


1.912 98(4) 


32 


924.7(16) 


1.299 30(7) 


1.528 98(9) 


-52.5(8) 


1.473 299(5) 


1.91018(2) 


64 


284.9(9) 


1.263 33(7) 


1.505 48(5) 


-11.0(4) 


1.473 543(2) 


1.909 374(9) 


128 


86.0(6) 


1.239 88(9) 


1.495 63(4) 


-2.45(19) 


1.473 5940(12) 


1.909 165(5) 


256 


25.5(3) 


1.22732(8) 


1.49199(2) 


-0.67(11) 


1.473 604 7(2) 


1.909 107(3) 


512 


6.9(2) 


1.222 04(6) 


1.490 859(16) 


-0.19(5) 


1.473 6077(4) 


1.909 090 7(15) 


1024 


2.11(11) 


1.220 24(4) 


1.490 538(9) 


-0.03(2) 


1.473 608 3(2) 


1.909 086 7(4) 


oo 





1.219 435. . . 


1.490 416... 





1.473 608 7... 


1.909 086 2. . . 



This simple analysis provides a practical way of approaching the thermodynamic limit 
without knowing the limiting position of the peak in advance. 

We first run a complete simulation for some small lattice, covering the whole range of 
rh. This provides a first approximation to the position of the peak. For growing lattices, 
we just compute two or three points at both sides of where we think the maximum is going 
to be. Our objective is not to reconstruct the whole peak of p{rh), just to find a good 
approximation to the point TOpe^k where {h)rh,(S vanishes. We use the same procedure 
as in Section 5.3, finding the zero of the cubic spline and interpolating the physical 
observables. Actually, if the position of the peak is sufficiently bounded we could just 
place one point very closely at either side and use a linear interpolation. 

With this procedure we are able to compute the tethered mean values of the relevant 
physical quantities at the peak with a minimum of numerical effort. Here we shall apply 
this method to the energy and we shall also characterise the approach of the peak to 
my. To the latter purpose, we have computed = (/i)m^,/3 for several lattice sizes and 
studied how fast it approaches zero. We also give the values for the position of the peak 
(Table 8). 

Following the above analysis we should find that 



^peak,/3 






L~ 




™peak,/3 


- my 


— 


■ L 






hy 


= Af, ■ 


L- 





with ( a 2. We have tried to fit these quantities to a power law for (3 = 0.4473, but 
we found its behaviour to be more complex. Instead, we present in Table 9 the result 
of taking the points for L and 2L and computing the effective exponent between them 
(that is, the value of ( so that the power law would pass exactly by those points). As we 
can see, our results are always C < 2, even though this exponent grows with L. 

We believe this was caused by the proximity of the critical point, so we ran analogous 
simulations for /3 = 0.6. We were able to complete this new computations in very little 
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Tabic 9 

Rate at which several obscrvablcs approach zero. Wc consider a functional form A ■ L^^ and compute 
the effective exponent C, from the ratio of the computed values at consecutive lattice sizes. Wc consider 
three exponents, C,f^, and C,u for the evolution of /ly, t'^i^^^-^ and ttp^j^j-, respectively, see equation (74). 
We observe that for (3 = 0.6 the effective exijoncnt aijpioaches 2, as expected from the discussion in the 
text, while for fi = 0.4473 the proximity of the critical point complicates the analysis. 







/3 = 0.4473 






(3 = 0.6 




L 










Crfi 


C« 


16 


1.690(3) 


0.6394(14) 


1.168(4) 


2.47(2) 


2.43(2) 


1.83(3) 


32 


1.699(5) 


0.864(3) 


1.356(6) 


2.25(5) 


2.24(5) 


1.93(5) 


64 


1.729(12) 


1.102(7) 


1.531(12) 


2.17(12) 


2.16(13) 


1.87(10) 


128 


1.75(2) 


1.375(16) 


1.73(3) 


1.9(3) 


1.89(13) 


1.9(2) 


256 


1.88(5) 


1.60(4) 


1.83(6) 


1.8(4) 


2.0(5) 


2.2(5) 


512 


1.71(9) 


1.70(8) 


1.85(12) 


2.7(12) 


1.4(10) 


3.1(12) 



time, following the above procedure. For example, for L = 512 the position of the peak 
was so well bounded that we just computed one point at either side. 

Comparing Table 8 with Table 2 we see that for = 0.6, with a computation effort 
almost 40 times smaller, we have obtained a result an order of magnitude more precise 
than what we had at /3c. Recomputing the effective exponents for these new simulations 
we obtain results compatible with C = 2. Notice that for this temperature the error in 
the exponents is much bigger than that for /3 = 0.4473. The reason is clear from Table 8. 
The left hand sides of equations (74) are now much closer to zero than in /3 = 0.4473, 
yet their errors are only slightly smaller. In the computation of the effective exponents 
only the relative errors matter, which explains our bigger uncertainties. Notice, however, 
that we have been able to distinguish values for /ly of order 10~^ from zero and that we 
have located the peak with seven significant figures. 

9. Conclusions and outlook 

We have presented the Tethered Monte Carlo method, a completely unspecific formal- 
ism to reconstruct the effective potential of the order parameter. This method is based 
on a new statistical ensemble, which we have described in detail. The tethered ensemble 
not only allows to reproduce the canonical results, but is also specially suited to the 
study of the broken symmetry phase and the effects of an external magnetic field. 

We have implemented this formalism in the D = 2 Ising model, were we can make 
exhaustive checks of our results, either against exact solutions or high precision compu- 
tations performed with canonical cluster algorithms. We have not have tried to optimise 
the efficiency for this particular model, choosing instead the most generally appliable 
method: a standard local Metropolis algorithm. The possibility of developing a cluster 
algorithm for the tethered formalism will be considered in the future. 

The method has been tested to a high degree of accuracy in a large variety of situ- 
ations: the critical point, the scaling paramagnetic region, the thermodynamic limit in 
the ferromagnetic region and the nonlinear response to an external magnetic fields. Our 
Monte Carlo implementation has proven to be remarkably efficient for the ferromagnetic 
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phase and specially for computations with an external magnetic field. 

One of its most conspicuous features is the absence of critical slowing down for all func- 
tions of the magnetisation, even though we have used a local algorithm. This is probably 
due to the fact that the tethered ensemble instantaneously propagates information to 
the whole lattice by means of a 'magnetic bath'. Our algorithm allows us to study each 
magnetisation independently without having to wander randomly in the magnetisation 
space, as is the case for multicanonical or Wang-Landau computations. 

We expect Tethered Monte Carlo to be of great help whenever large tunnelling barriers 
appear, associated to observablcs not present in the Hamiltonian. A non exhaustive list 
of instances are the standard magnetisation in the Random Field Ising Model [44], the 
staggered magnetisation for the Diluted Antiferromagnet in a Field [45] or the Polyakov 
loop for lattice gauge theories [8]. In these cases 'tethering' extensive quantities other 
than the magnetisation could prove an invaluable guide for the exploration of phase 
space. In particular, we wish to apply this formalism in the near future to undertake a 
thorough study of the condensation transition [21,46]. 
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Appendix A. Some numerical details 

A numerical implementation of equation (28) can be done in several different ways of 
equivalent numerical accuracy. Here we briefly explain our choices. 

Recall that we have represented {h)m,f3 with a cubic spline interpolation. We have 
not used the so called natural spline, which imposes vanishing curvature for {h)m,j3 at 
mmax and mmin- Instead, we have estimated the derivative of this function at both 
ending points with a parabolic interpolation. To compute the canonical average of (28) 
we also represent {0)m,f3 with a cubic spline. However, naively applying this interpolation 
scheme for the exponential, p(m) = ex.p[NnN{m)], could introduce strong integration 
errors. Fortunately, this can be easily solved by accurately representing the exponent 
Nf2]\[, which is a smooth function (recall the curve in logarithmic scale of Fig. 1). 

i7jv(™, /?), being the integral of the cubic spline for {h)rh.p, is a fourth order polynomial 
between each pair of consecutive points in the m grid, which can be exactly computed 
from the coeflacients of the cubic spline. To avoid losing precision, we evaluate J7jv("i,/3) 
at an extended grid that includes 3 equally spaced intermediate points between each pair 
of simulated values of m. This way the Lagrange interpolating polynomial for c;ac;li seg- 
ment of the extended lattice (two original points plus three intermediate ones) represents 
the exact integral of our spline. 

Of coiuse. the pdf, cxp[N f2j\[{m, jS)], is not a polynomial anymore. It is, however, a 
smooth function so a self consistent Romberg method [31] provides an estimate of the 
integral (28) with any required numerical accuracy. Notice that this yields the basically 
exact results for a given interpolation of {h)rh,0, but it does not cure any discretisation 
errors introduced by the spline. 
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Typically, even with a very moderate effort, the Romberg integration error has been 
much smaller than the statistical one. There is one exception: the fluctuation-dissipation 
formulas, such as (7), because of the large cancellations between the two terms. To solve 
this problem, we have computed the fluctuation-dissipation formula as a sum of two 
squares: 



N-^C = {u')0 - {u)l = y dm p{m) [{u'U ' + ' (") 



dm p{m) ((u - {u)rhf)rh + {{u)m - {u)py 



(A.l) 



In spite of this, as a consistency check, we have also employed the original equation and 
forced the Romberg integral to yield the same value, by reducing its tolerance. 
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